講義資料 計算生命科学の基礎Ⅳ （→終了しました） | 計算科学教育センター

QM/MM法を用いた

タンパク質の機能解析

広島市立大学大学院情報科学研究科

鷹野 優

遠隔インタラクティブ講義

「計算生命科学の基礎（IV）」

2017年11月22日

概要

I.  序論

II.  QM/MMの理論

III. 金属タンパク質中の補欠分子の電子構造

III-1. 酸素運搬タンパク質ヘムエリスリン III-2. 光合成反応中心のスペシャルペア

IV. 生体高分子生化学的機能解析のための分子計算技術の開発

IV-1. 半経験的分子軌道法によるQM部分の高速化

IV-2. 超並列QM/MM-MD連成計算プログラムPlatypus

タンパク質は生命活動を担う

タンパク質 ＝ 生命活動をおこなう分子

免疫

生きるのに必要な様々なはたらきをする

不具合が起きると病気に！

Molecule of the Month © David S. Goodsell and RCSB PDB licensed under CC-BY-3.0 license

タンパク質と病気

通常骨折にならないようなごく軽微の外的要因でも骨折する。

ヒトの組織の構造タンパク質、I型コラーゲンの変異が原因

Gly

Ala

正常

骨形成

全症

Gly

_→

Ala

構造を維持できなくなる→骨形成不全症

タンパ

質

一部が変化

https://sites.google.com/a/comp-prot-sci.org/tanpakushitsu-no-rittai-kouzou-nyuumon/chapter-6 Molecule of the Month © David S. Goodsell and RCSB PDB licensed under CC-BY-3.0 license

タンパク質と薬

抗インフルエンザ薬：タミフル

Our goal

分子構造

機能

構造機能相関

分子のかたち

電荷分布 軌道相互作用

（例）Croタンパク質

DNAと強く結合し、RNAポリメラーゼによる転写を邪魔します DNA

Croタンパク質 赤：負電荷 青：正電荷

Our goal

分子構造・電子構造は機能発現を

どのように制御しているか？

電子構造

分子構造

機能

分子構造ー電子構造ー機能相関

分子のかたち

タンパク質の電子構造

1) タンパク質は巨大で、ヘテロで、ダイナミックな系

2) タンパク質の機能の発揮は局所部位の電子状態に依存する

光合成反応中心

タンパク質のつくる環境は電子構造に影響

タンパク質の電子構造

1) タンパク質は巨大で、ヘテロで、ダイナミックな系

2) タンパク質の機能の発揮は局所部位の電子状態に依存する

タンパク質のうごき(ダイナミクス)は電子構造に影響

1) タンパク質は巨大で、ヘテロで、ダイナミックな系

2) タンパク質の機能の発揮は局所部位の電子状態に依存する

活性中心そのもの

活性中心をとりかこむタンパク質 溶質分子をとりかこむ溶媒分子

アンサンブル ダイナミクス

計算サイズ 小

大

資源（計算機・人間の能力）には限りがある

「どのようにしてモデル化するか?」

タンパク質の電子構造

功：計算量の減少

  系の単純化

罪：情報の欠損

モデル化の功罪

資源の分配をどうするか

なにを理解したいのか？

タンパク質の電子構造

QM/MM法の適用

QM

MM

Electrosta+c

effect

Steric

effect

Photoactive Yellow Protein

1) タンパク質は巨大で、ヘテロで、ダイナミックな系

QM/MM法

h0p://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/popular-chemistryprize2013.pdf

・多階層連結計算法の一つ

・電子構造が重要な部分（化学反応）は量子力学(QM)を

 それ以外は分子力場(MM)を用いる

・1976年にWarshelとLevittが提案

QM/MM法

Inner

(QM)

Outer

(MM)

_System

System = Inner + Outer + (Link)

共有結合

・全系(System)をQM領域(Inner)とMM領域(Outer)に分割

・QMとMMの境界に共有結合が含まれる場合は

QM/MM法に使う計算手法

・基本的に制限はない

QM: 密度汎関数法(DFT)

   Hartree-Fock法(HF)

   post-HF法(MP2, CASSCFなど）

  半経験的分子軌道法（AM1, PM6法など）

QM/MM法の種類

EQM/MM sub

System

(

)

(

)

(

)

(

)

全系のエネルギーの表現法に関して大きく二種類ある

減法表現（ONIOM法）

加法表現（QM/MM法）

EQM/MM add

System

(

)

(

)

(

)

(

)

+

‒

=

=

+

+

相互作用 QM-MM

タンパク場の効果

mechanical embedding

electrostatic embedding

QMの電子密度はMMと相互作用しない

Steric effectのみ考慮

QMのハミルトニアンにMMとの静電相互作用をとりこむ

Electrostatic effectとSteric effectを考慮

MMの電荷がQM領域に近いときには過分極を起こす →・電荷スケーリング法

   QM境界付近のMM電荷の大きさを調整  ・電荷シフト法

QM/MM境界問題

QM/MM境界に化学結合が存在する場合

QM空間の切断->不対電子対の発生->QM（分子軌道）の歪み

正しいエネルギー及び力が得られない

1.

どのように切るか？

QM/MM境界問題

1. どのように切るか？

(1)Link atom method

(2)局在化軌道法

QM

MM

Hydrogen Atom etc…

不対電子対

QM領域とMM領域の境界にある 化学結合の末端を水素原子などで 置き換える

QM領域とMM領域の境界にある 化学結合を局在化軌道(SLBOs)を 用いて表現する。

QM/MM境界問題

(2) 線形応答関数解析

δρ

(

_r

)

=

δρ

(

r

)

δ

v

₍

_r

!

₎

δ

v

₍

_r

!

₎

∫

d

_r

_!

QM/MMモデル化に起因する

QM領域の電子構造の誤差の最小化

2. どこを切るか？

(1) 経験的、化学的洞察

線形応答関数解析

QM/MMモデル由来のQM領域の誤差  の最小化

従来：QM/MMモデルによる誤差 の最小化

δ

v

₍

_r

!

₎

δρ

(

_r

)

(i)  共有結合では電子トランスファーが大きい平衡核間距離ではなく、

結合不安定が起こる(結合が切れかけた)状態で最大となる。

(ii)  トランスファーブロックがない状況ではサイト数が多くなる

につれフリーデル振動型の振動・減衰となる。

(iii) アラニンジペプチド, グルタチオン, アラニントリペプチド系 ではsp3_{接合点が伝播をブロックする。水素結合経由の伝播} は約0.001(ex. 0.1 Hartreeの に対し0.0001の誤差伝播)。

  αヘリックス系では明白に誤差伝播が構造化される(右図)。 δρiatom δvjatom =1 δρ_iatom δv_jatom =−₁

δρ_iatom δv_jatom =0.0

RGB=(1-LRF, LRF, 0.5)/max LRF

δv_(r!)

[ガイドライン]

Ueda et al. Int. J. Quantum Chem. 113, 336–341 (2013).

QM/MM境界問題の定量化が可能に！

QM/MM法による構造最適化

計算コスト

 QM > QM/MM >> MM

(1)断熱アプローチ

QMのmacroステップではMM環境は完全に   構造最適化される

(2)交互スキーム

QM部分とMM部分の構造最適化に異なるアルゴリズムや 座標系を使う

 QM: macroiteration（内部座標＋擬ニュートンアルゴリズム）  MM: microiteration（カーテシアン座標＋共役勾配法）

microiterative scheme

溶媒分子のとりあつかい

多数の溶媒分子の配置・配向をどのように表現するか？

      (112アミノ酸残基） 水 24,039原子（8,013分子）

系のほとんどが水分子！

水分子の配置＋配向

（高い自由度）

①

②

構造の異方性

水素結合

溶媒分子のとりあつかい

多数の溶媒分子の配置・配向をどのように表現するか？

①

②

構造の異方性

水素結合

1. あらわにとりあつかう

・古典力学 MM

2. モデル化してとりあつかう

・統計モデル

RISM-SCF, ER

・連続誘電体モデル

PCM, CPCM, SMx, COSMO

古典力学 (QM/MM, ONIOM)

C O

O CH₃

CH_{3 H} O _H

O H H O H H O H H O H H O H H O H H O H H O H H C O O CH₃

CH_{3 H} O _H

O H H O H H O H H O H H O H H O H H O H H O H H MM

活性中心：量子力学、 タンパク質：古典力学

溶媒環境：古典力学

利点：精度

統計モデル (RISM-SCF, ER)

C O

O CH₃

CH_{3 H} O _H

O H H O H H O H H O H H O H H O H H O H H O H H

溶媒の分布関数

活性中心：量子力学、 タンパク質：古典力学

溶媒環境：統計理論

利点：精度・計算量

欠点：ダイナミクスのとりこみ

連続誘電体モデル

(

PCM

, CPCM, SMx, COSMO)

C O

O CH₃

CH_{3 H} O _H

O H H O H H O H H O H H O H H O H H O H H O H H C O O CH₃ CH₃ Cavity water ε = 78.39

ε = 1

活性中心：量子力学、 タンパク質：古典力学

溶媒環境：連続誘電体

利点：計算量

欠点：精度（水素結合のとりこみ）

PCMの形式

H

=

H

0

+

V

ハミルトニアン

真空中の分子（溶質）の ハミルトニアン

溶質−溶媒相互作用の ポテンシャル

溶質分子が溶媒分子に囲まれたとき、

PCMでは、

1. 溶媒分子を誘電率εの連続媒質におきかえ、

PCMの形式

H

=

H

0

+

V

ハミルトニアン

真空中の分子（溶質）の ハミルトニアン

溶質−溶媒相互作用の ポテンシャル

溶質分子が溶媒分子に囲まれたとき、

PCMでは、

2. 溶質の形にもとづいてcavity（空孔）をつくり、

ε = 78（水）

PCMの形式

H

=

H

0

+

V

ハミルトニアン

真空中の分子（溶質）の ハミルトニアン

溶質−溶媒相互作用の ポテンシャル

溶質分子が溶媒分子に囲まれたとき、

PCMでは、

3. cavity（空孔）の中に溶質分子をおく

ε = 78（水）

ε₀ = 1.0

H O H H O H H O H

ε₀ = 1.0

PCMの形式

H

=

H

0

+

V

V

int

=

∫

ρ

( )

r

Φ

( )

r

d

r

真空中の分子（溶質）の ハミルトニアン

溶質−溶媒相互作用の ポテンシャル

溶質の電荷密度 溶質によって分極した溶媒（連続誘電体）の つくるcavity内の電場

溶媒効果を溶質の電荷密度と溶質によって分極した連続誘電体の つくる電場との相互作用として表す(Reaction field法(RF法))

H O H H O H H O H

ε₀ = 1.0

ε = 78（水）

PCMの形式

面積素片k（位置_rk,電荷Q_k, 面積S_k)

Φ

( )

_r

=

Q

r

−

r

∑

誘電体のつくる電場は、cavity表面近傍に分布する分極電荷の与える電場に等しい 計算では、cavity表面を面積素片(tessera)にわけて、

面積素片の位置ベクトル 面積素片のもつ電荷

nk

面積素片kの 法線ベクトル

Φ_in

( )

Φout

( )

Q

=

ε

−

1

4

πε

S

k

(

∇Φ

( )

r

•

n

)

面積素片のもつ電荷は古典電磁気学から求められる

比誘電率 面積素片の面積 _{面積素片の法線ベクトル}

PCMの形式

Φ

( )

_r

=

Φ

( )

_r

+

Φ

( )

_r

溶質の電荷密度により つくられるcavity内の電場

溶質によって分極した溶媒（連続誘電体）の つくるcavity内の電場

cavity内の電場は溶質によるものと溶媒によるものにわけることができる

Φ_in

( )

Q_k

Φ_in,solvent

( )

Q_{k の計算には が必要}

Φ_in

( )

( )

くりかえし計算で解く

生物と金属

多量

Na

K

Mg

Ca

微量

Fe

Zn

Cu

Mn

V

超微量

Ni

Cr

Co

Mo

W

バランスを保つことで生物の生命維持に必須のはたらきをする

呼吸・光合成・窒素固定

金属タンパク質

金属をもつ生体分子の特徴

（環境にやさしい触媒）

・ハーバーボッシュ法

例：窒素（反応しにくい安定な気体）からアンモニアを生成

N

+ 3H

→ 2NH

N

+ 3H

2NH

高温・高圧

(500 ℃, 300気圧)

・ニトロゲナーゼ（鉄：Fe, モリブデン：Mo)

N

+ 3H

2NH

常温・常圧

金属補欠分子

ヘム (Fe)

スペシャルペア (Mg) [4Fe-4S]クラスター (Fe)

酸素発生複合体(Mn)

タンパク質の生物活性において重要な、タンパク質に結合する非タンパク質

Inside-out アプローチ

（例）シトクロムc酸化酵素のヘム a in CcO

1. 金属補欠分子固有の電子構造を調べる

Takano and Nakamura, J. Comput. Chem. (2010)

2. タンパク質中での電子構造を調べる

static effects (構造歪み、静電相互作用)

dynamic effects (ゆらぎ)

3. 1と2を比較する

タンパ 質 電荷移動 を強める

タンパク質なし タンパク質あり

PropA 電荷変化が

遷移金属と生体分子の相互作用

電荷

軌道

スピン

強相関電子系

酸化数の変化

スピン状態の変化 stepwiseな電子移動

軌道・電荷・スピンが協同的に影響しあう

強相関電子系

金属イオンの電子間にはたらく強いクーロン反発相互作用

さまざまな電子配置・電子相関を考慮にいれないといけない （例）酸素発生複合体 (Mn4CaO5クラスター)

可能な電子配置数 Mnの3d軌道

Oの2p軌道 約100京通り！

ab initio

DMRG法で可能になった

Density Functional Theory (DFT)

In DFT

E

Ts : kinetic energy of non-interacting system

E

XC : exchange-correlation energy

E

₌

E

T

+

E

+

E

+

E

Approximation

E

DFT

=

E

+

E

+

E

+

E

+

E

Exchange-Correlation Energy

DFTの精度は交換相関汎関数項に依存する

E

= (

E

–

E

) + (

E

–

E

)

Correction of kinetic energy

Correction of electron-electron interaction

Density Functional Theory (DFT) is one of the most popular tools

performance ex. Mean absolute error of heat of formation for the G2 dataset (kcal/mol)

DFT 3.0 (B3LYP), 7.1(BLYP), 90.8 (SVWN) HF 149.2

MP2 23.8 CC 11.5

cost ex. Benzene (C₆H₆) single point

DFT O(N3_)~O(N4_{) 1 min 38 sec (SVWN), 2 min 20 sec (BLYP)}

HF O(N4_{) 1 min 48 sec}

MP2 O(N5_{) 2 min 36 sec}

CC O(N7_{) 2 hours 29 min 20 sec}

Xu, X. and Goddard III, W. A. , Proc. Natl. Acad. Sci. U.S.A. 2004_{, 101, 2673-2677.}

Jensen, F. Introduction to Computational Chemistry, 2nd ed. Wiley (2007)

Excellent performance-to-cost ratio !

Obscure systematic improvability

There is not clear path to the exact solution of the Schrödinger eq., since DFT show very good accuracy for a wide range of applications due to the error cancellation.

DFTの欠点

Obscure systematic improvability

There is not clear path to the exact solution of the Schrödinger eq., since DFT show very good accuracy for a wide range of applications due to the error cancellation.

III-1.酸素運搬タンパク質ヘムエリスリン

ヘムエリスリン

Hr is an oxygen transport protein found in several marine invertebrates.

O₂

DeoxyHr OxyHr

FeII O FeII OO His77 His73 His101 His54 His25 O O Asp106 H Glu58 FeIII O FeIII OO His77 His73 His101 His54 His25 O O Asp106 Glu58 O O H hydroperoxide J

ab = –13 cm–1 Jab = –77 cm–1

[Fe(µ-OH)Fe] [Fe(µ-O)Fe]

In the reac+on with molecular oxygen, two electrons from the FeII_–FeII _{center and a}

二核鉄活性中心

Diiron ac+ve sites with µ-O and/or µ-OH bridges are found in many proteins.

The structure of the ac?ve site are similar to each other,

but the func?ons are different.

Fe O Fe OO His His His His O O Asp His Glu O O H Fe O Fe O Glu His Glu His Glu O Glu O H H

H₂O

Oxygen transport Oxygen ac?va?on

(Conversion of methane to methanol)

Metane monooxygenase (MMO)

Hemerythrin (Hr)

Fe O Fe

O O His His Glu O Glu Glu H₂O

H₂O O

Asp

Ribonucleo?de reductase (RNR)

Oxygen ac?va?on

以前の研究

E

α β

Opt Fe2

d_xy(δ_⊥)

dx2_–y2(δ)

dxy(δ_⊥) dyz(π) d_zx(π_⊥) dz2(σ)

α β

α β

X-ray Fe2

dx2_–y2(δ)

dxy(δ_⊥) dyz(π)

dzx(π_⊥) dz2(σ)

dxy(δ_⊥)

O₂ –8.0 –0.8 –9.0 –10.0 –11.0 –12.0 –13.0 –14.0 –15.0 –16.0 π∗,

π∗,π_⊥∗

π_⊥∗

Takano, Y.; Isobe, H.; Yamaguchi, K. Bull. Chem. Soc. Jpn. 2008, 81, 91–102

タンパ 質環境 活性中心 構造を歪めるこ

目的

ヘムエリスリンのタンパク質環境による

静電相互作用は

アイデア

タンパク質環境

タンパク質環境を構成

する原子電荷との

直接の静電相互作用

活性中心の電子状態に

よって誘起される

原子電荷の揺らぎ

（分極効果）

＋

＝

手法の選択

How do we determine an appropriate broken-symmetry theoretical treatment?_"

1. Use a method which many people use. → UB3LYP

UB3LYP provides very accurate results for organic systems. But how about transi?on metal systems?

2. Use a method which can reproduce experimental data.

Magne?c coupling constants (J

ab) were reported.

(OxyHr: –77 cm–1_{, DeoxyHr: –13 cm}–1₎

化学結合と磁気的相互作用

Magnetic Coupling Constants,

J

_ab

a

H

₌

₋

₂

J

_ab

_{S S}

b

⋅

∑

Heisenberg Hamiltonian_"

Fe O Fe OO O O a

S

S

Magnetic coupling constants, J

ab, can be calculated using the Heisenberg hamiltonian. " Each iron site is regarded as a spin site._"

◆ J_ab values!

J

ab

=

LS

E

₋

E

HS

S

₋

LS

S

ab > 0: ferromagnetic"

J

計算方法

◆ Electronic Energy"

Kinetic energy term_"

Nuclear-electron attraction_" Coulomb integral_"

Exchange-correlation functional_"

E ₌ ET _{+ E}V _{+ E}J _{+ E}XC

! ET ! EV ! EJ !

EXC_!

EXC

= C

1EHF + C2ESlater + C3EBecke + C4EVWN + C5ELYP !

C

1 C2 C3 C4 C5

UHF 1.0 0 0 0 0

UB2LYP 0.5 0.5 0.5 1.0 1.0

UB3LYP 0.2 0.8 0.72 1.0 0.81

UB3LYP* 0.15 0.85 0.72 1.0 0.81

モデル

DeoxyHr (PDB ID: 1HMD)_"

OxyHr (PDB ID: 1HMO)_"

His_"

Asp, Glu_"

Magnetic Coupling Constants,

J

_ab

Solvation UHF UBHandHLYP UB3LYP UBLYP

真空 –3.49 –10.4 –145.1 –396.1

PCM –3.59 –10.6 –128.2 –326.4

PCM + MM –3.58 –10.6 –137.9 –398.8

DeoxyHr (Exptl.: –14 cm

₎

Magnetic Coupling Constants,

J

_ab

Solvation UHF UBHandHLYP UB3LYP UBLYP

真空 –14.1 –49.0 –145.1 –396.1

PCM –17.2 –64.0 –128.2 –326.4

PCM + MM –16.7 –64.4 –137.9 –398.8

OxyHr (Exptl.: –77 cm

₎

Orbital Energy Gap

タンパ 質環境 よる静電相互作用 活性中心

Orbital Energy Gap

O

π*

軌道

‒0.8

‒8.84

E

(eV)

Fe

dδ

軌道

‒8.27

‒6.71

Opt X-ray X-ray+PCM_+MM

タンパク質による構造歪み × ○ ○

タンパク質による静電相互作用 × × ○

まとめ

BHandHLYP法はOxyHrとDeoxyHrの双方で妥当的な 磁気的相互作用を与えた

ヘムエリスリンの蛋白質環境による静電相互作用は 活性中心の金属のd軌道のエネルギー準位を制御し 酸素結合に寄与している

1.

III-2. 光合成反応中心スペシャルペア

光合成

植物などの光合成色素をもつ生物がおこなう、

光エネルギーを化学エネルギーに変換する生化学反応

ATP

+

h

ν

二酸化

炭素

(CO

)

水

(H

O)

酸素

(O

)

糖

(C

H

O)

+

+

光合成電子伝達反応

光合成電子伝達反応（明反応）

Essen+al細胞生物学

光化学系

反応中心: 光エネルギーを電子エネルギーに変換

光合成反応中心

chromophores

Special pair (SP)

Bchl_L

Bchl_M

Pheo_L Pheo_M

Q_L Q_M

Fe

h

ν

e–

光合成反応中心

His(L173)

SP_L

SP_M

His(M202)

Special pair (SP)

L-branch

M-branch

Bchl_L Bchl_M

Pheo_L Pheo_M

Q_L Q_M

Fe

h

ν

色素が疑似二回対称性を持つにもかかわらず 電子移動経路は一方向

スペシャルペアカチオンラジカルの

スピン密度

Lubitz et al., Acc. Chem. Res. 2002, 35, 313

スペシャルペアカチオンラジカルのスピン密度は非対称

Leu(M160)His Leu(M160)His + Phe(M197)His Arg(L135)Glu Arg(L135)Leu Phe(M197)His Arg(M164)Leu Arg(M164)Glu Leu(L131)His His(L168)Phe wild type

SP_M SP_L Leu(M160)His + His(L168)Phe Leu(M160)His + Leu(L131)His

His(L168)Phe + Phe(M197)His Leu(L131)His + Phe(M197)His Leu(L131)His + His(L168)Phe

100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0%

SP_L

SP_M

目的

スペシャルカチオンラジカルの

スピン密度の非対称性の原因は何か？

タンパク場の効果

l  スピン密度の偏りが再現された

l  スペシャルペア単体でもスピン密度に偏りが生じている

l  タンパク質の静電相互作用がスピン密度の非対称性が強める

Yamasaki, H.; Nakamura, H.; Takano, Y. Chem. Phys. Lett. 2007, 447, 324–329

[QM]方法:UB3LYP 基底関数6-31G(d)(0.25) [MM] AMBER96

側鎖の影響

Acetyl group

13keto carbonyl group

Nonpolar groups

Methyl ester (Mes) group

Phytyl-II (Phyt) group

Phytyl-I (Phyt) group

Mes基とPhyt基の配向の違いにより

側鎖の影響

Yellow: positive Δρ

Blue: negative Δρ

Mes

基

Phyt

基

配向が異

る

S244(L)

Gly280(M)

Cys247(L)

Gly283(M)

配向の違いは保存されている

1aig 2i5n 1eys 2axt

Rhodobacter sphaeroides

Rhodopseudomonas viridis

Themochro- matium tepidum

Thermosynecho- coccus elongatus

タンパク質環境が配向の違いを与える

Side view

M-side L-side

Top view Phe181

Leu131 Ala176 Ser244 Met248 Tyr210 Leu209 Ser206 Gly280 Ile284 Leu156 Thr277

Rb. sphaeroides (red)

Rh. viridis (blue) T. tepidum (green)

T. elongatus (yellow)

Phe180

４つの結晶構造の重ね合わせでは、

配向の違いを与える

アミノ酸残基は保存されている

Ile/Ile

Ser/Thr Leu/Val

Ala/Gly

配向の違いを与える

アミノ酸残基は保存されている

?/Leu

Leu or Phe/Val Leu/Val Tyr/Leu

Gly or Ala/Leu

Ile/Ile

Ser/Gly

NMRによる検証

まとめ

1. Mes基とPhyt基の配向の違いがスペシャルペアの   電子構造の非対称性を引き起こした

2. タンパク質の静電相互作用がスピン密度の非対称性を強めた

Inside-out アプローチ

1. 金属補欠分子固有の電子構造を調べる (QM) 2. タンパク質中での電子構造を調べる (QM/MM)

static effects (構造歪み、静電相互作用)

dynamic effects (ゆらぎ)

3. 1と2を比較する

高精度化した半経験的分子軌道法による

QM部分の高速化

+

超並列

ab initio

QM/MM-MD連成計算プログラム

IV-1. 高精度化した

半経験的分子軌道法による

QM部分の高速化

Saito, T.; Kitagawa, Y.; Takano, Y. J. Phys. Chem. A 2016_{, 120, 8750–8760.}

size of QM region UDFT/DZP level

600

number of steps in MD run

e.g.,) 1ns = 1_×106 _step

it adds up to the 1,000,000 QM/MM computations time for QM/MM

single point calculation / s

Provided that every QM/MM energy and force computation takes 10 min.,

→ one MD run takes 20 years !!

QM/MM-MD

0.1

1ns MD run requires 70 days. number of steps in MD run

e.g.,) 1ns = 1_×106 _step

Semiempirical MO (SE-MO)

roughly 2 or 3 orders of magnitude faster

size of QM region UDFT/DZP level

600

M. J. S. Dewar and W. Thiel, J. Am. Chem. Soc. 99, 2338 (1977). W. Thiel, WIREs. Comput. Mol. Sci. 4, 145 (2014).

time for QM/MM

single point calculation / s

L. Goerigk and S. Grimme, J Chem. Theory Comput. 6, 107 (2010). L. Goerigk and S. Grimme, Phys. Chem. Chem. Phys. 13, 6670 (2011).

P. Dral et al., J. Chem. Theory Comput. 12, 1097 (2016).

Given the applicability to a wide range of molecules,

we chose PM6, which supports 70 elements in the periodic table.

MAEs (kcal/mol)* 遷移金属への対応

MNDO 27.3 No

AM1 14.7 No

PM3 11.3 No

PM6 10.2 YES

OM3 6.3 No

GMTKN30-CHNOF database

(430 reactions mostly involving closed-shell molecules)

J. J. P. Stewart, J. Mol. Model. 13, 1173 (2007).

＊_{mean absolute errors,excluding the MB08-165 subset}

Some problems using SE-MO?

■ cyclobutadiene (1C₄H₄) radical

non-planar linear with the H-C-H angle of 180˚

■ transition state of C=C epoxidation with a Mn complex

fail to find a TS

×

C

H H

C

C C

C H H

H H

The spin-unrestricted PM6 (UPM6) computation does not work.

bent : H-C-H = 102.3˚

目的

the original PM6

reparameterized PM6 (rPM6)

The number of optimizable electronic and core−core parameters adds up to

43 (P) for the basic elements (H, C, N, O).

••••

Our primary purpose is to improve the accuracy of PM6 in describing electronic and geometric structures of open-shell species,

for which a slight modification was made.

training set containing (di)radicals using

the UB3LYP/6-31G* model as reference data

part of GMTKN30 was also included to maintain the performance for general-purpose applications

Application to Methylene

(included in the training set)

Method 1_CH

2 3CH2 S-T gap

C-H H-C-H C-H H-C-H [kcal/mol]a

UAM1 1.068 140.6 1.062 151.7 34.1

UPM3 1.067 131.3 1.063 149.9 42.7

UPM6 1.018 180.0 1.019 180.0 31.5

UrPM6 1.081 111.9 1.039 152.6 11.2

UB3LYP/6-31G* 1.090 113.5 1.082 133.1 6.0

Expt. 1.107 102.3 1.077 134.0 9.0

a _{With the AP correction}

It is not surprising that the present UrPM6 shows great improvement over the original UPM6.

rPM6 shows great improvement

over the original UPM6

■ transition state of C=C epoxidation with a Mn complex ■methylene (1,3 CH2) radical

■ cyclobutadiene (1C4H4) radical C

H H

C

C C

C H H

H H

bent

H-C-H = 111.9˚

square planar (D_4h) planar

CPU time ratio (vs. DFT) _"

force calc. 1 : 4608_"

Hessian calc. 1 : 3559_"

IV-2. 超並列

ab initio

QM/MM-MD

連成計算プログラム

Platypus

超並列

ab initio

QM/MM-MD連成計算プログラム

Platypusの開発

(PLATform for dYnamics Protein Unified Simulation)

E

=

_E

+

_E

+

_E

+

Chain-of-State法 (疎結合並列計算）

Platypus-QM Platypus-MM

並列性能に優れた新規

長距離ポテンシャルの開発 (zero-multipole法)

HF, DFT, CASSCF, CIS, MP2プログラムの

現在までの研究開発成果

• RDF, UHF, R-DFT, U-DFT, CIS, CIS(D), CIS-DFT, MP2, CASSCF の実装

• エネルギー・力計算における積分計算のハイブリッド並列化とSIMD化

• 2電子積分計算の原子軌道基底から分子軌道基底への変換部分の高度化

• CASSCF計算におけるdirect CI 計算の高速化（CAS(16,16)まで可能）

QM部分の高速化・超並列化

・ 粗結合 (chain-of-state法)を利用した超並列化

「京」での並列性能(QM)

並列化率 99.9888%

光合成反応中心 スペシャルペア

(280原子）

RHF/cc-pVDZ

 総数32,768コア

 並列化率 99.9888%

 実行効率 7.27∼2.24％  (8192コア：4.26%)

CASCI(16,16)/6-31G**

 総数16,384コア

 並列化率 99.9728%

 実行効率 15.43∼6.80％  (8192コア：7.85%)

!" # !" #$ #$# $ %

「京」での並列性能(QM/MM-MD)

98,304コアのハイブリッド並列の実行

Ace-ALA-NMe +

1,080 水分子

Strong scaling

 総数2,048コア

 並列化率： 99.7893%

 実行効率 10.60∼1.94％ RHF/cc-pVDZ

Weak scaling (Chain of State)

 総数98,304コア

 128プロセス x 768レプリカ  実行効率 3.42∼3.20％

 Weak scaling効率：93.36%

励起状態QM/MDシミュレーション

sirius発色団(1,374原子)

SA-CASSCF(8,8)/MINI-4 -643.5 -643.5 -643.4 -643.4 -643.4 -643.3 -643.3 -643.2 -643.2

0 5 10 15 20 25 30

Ground state

1st excited state

T o ta l e n e rg y / a .u .

Time / ps

QM: 発色団 MM: 水

sirius発色団

タンパク質のQM/MM計算

1) タンパク質は巨大で、ヘテロで、ダイナミックな系

2) タンパク質の機能の発揮は局所部位の電子状態に依存する

「どのようにしてモデル化するか?」

功：計算量の減少 vs. 罪：情報の欠損

高精度化した半経験的分子軌道法によるQM部分の高速化

超並列ab initio QM/MM-MD連成計算プログラムPlatypusの開発

・QM/MM法の適用

Protein as amplifier

・光合成反応中心のスペシャルペア

タンパク質が電子非対称性を強める

タンパク質なし タンパク質あり

実験値

タンパク質環境の電気的・立体的摂動によって、補欠分子の 固有に持つ性質が強められ、制御されている

・ヘムエリスリン

タンパク質が

FeとO₂の相互作用を

タンパク質

包丁

料理人

活性中心

h0p://www.michiba-shunsara.jp/menu/sz044.html

溶媒環境

ヘムエリスリン 光合成反応中心

・中村 春木 （大阪大） ・山崎 秀樹 （大阪大） rPM6開発

・齋藤徹 （広市大） ・北河康隆 （大阪大）

共同研究者

Platypus開発

・中村 春木 （大阪大） ・中田一人 （NEC） ・米澤康滋 （近畿大） ・山中秀介 （大阪大）

外部資金